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Abstract 


Reliable  methods  for  accurately  predicting  neutron  production  cross  sections  are 
needed  for  a  variety  of  applications  including  space  radiation  protection  and  the  use  of 
heavy  ion  beams  in  physics  and  medicine.  A  three-step  model  is  utilized,  including  a 
modified  optical  model  for  abrasion  and  ablation  and  a  coalescence  model  for  neutron 
removal  in  the  production  of  deuterons,  tritons  and  other  heavier  nuclei.  Estimates  of 
spectral  and  angular  distributions  of  neutrons  are  in  good  agreement  with  published  cross 
section  measurements. 
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Chapter  1 


Introduction 


Introduction 

The  collision  of  atomic  nuclei  traveling  at  or  near  relativistic  speeds  with 
stationary  target  nuclei  produces  a  shower  of  nuclear  fragments  in  all  directions  and  at  a 
wide  range  of  energies,  including  energies  significantly  higher  than  those  of  the  original 
projectile.  Sources  of  these  high-energy  projectile  nuclei  include  space  radiation  of 
cosmic  or  solar  origin,  accelerators,  and  radiation  beams  used  in  medical  treatment.  Of 
particular  concern  is  the  production  of  secondary  neutrons  in  bulk  matter,  such  as  the 
shielding  and  components  of  manned  spacecraft  or  space-based  facilities  and  overlying 
body  tissue  that  shields  underlying  organs.  Some  estimates  of  dose  equivalent  for  Lunar 
or  Martian  base  scenarios  indicate  that  up  to  50%  of  the  dose  equivalent  behind  thick 
shields  may  come  from  secondary  neutrons  (Simonsen,  1991).  They  are  also  a  major 
source  of  radiation  exposure  for  flight  crews  in  high-altitude  aircraft  (NCRP,  1995). 
Because  it  is  not  possible  to  determine  experimentally  the  neutron  production  spectra  for 
every  possible  combination  of  energy,  projectile  and  target,  it  is  desirable  to  produce  a 
validated  model  that  can  be  used  to  provide  cross  section  predictions  as  needed  for  these 
applications. 

Objective  of  This  Work 

The  objective  of  this  work  is  to  further  develop  existing  theories  for  nucleon 
production  in  high-energy  nucleus-nucleus  collisions,  which  accurately  predict  neutron 


production  cross  sections  in  these  reactions.  The  expanded  model  will  be  based  on  the 
abrasion-ablation  model  of  Cucinotta,  Wilson,  and  Townsend  (Cucinotta,  1997)  as 
extended  to  include  a  coalescence  model  (Nagamiya,  1981),  which  accounts  for  the 
formation  of  composite  particles,  such  as  deuterons,  in  the  debris  of  these  high  energy 
collisions.  This  model  is  used  to  produce  spectral  and  angular  cross  sections,  which  are 
compared  to  published  experimental  data. 

Organization  and  Scope  of  Work 

Chapter  2  briefly  explains  the  history  and  background  of  heavy  ion  collisions,  as 
well  as  the  desirability  of  an  accurate  model  to  predict  neutron  production  in  these 
collisions. 

Chapter  3  provides  details  of  the  theory  behind  each  stage  of  the  model,  including 
modifications  and  assumptions  used  to  implement  the  model. 

Chapter  4  presents  the  results  of  the  model  and  detailed  comparisons  with  prior 
work  and  with  experimental  data. 

Chapter  5  contains  a  summary  of  this  work,  a  discussion  of  the  advantages  and 
limitations  of  the  model,  and  a  detailed  description  of  work  that  should  be  considered  in 
the  future. 

The  appendices  include  a  more  rigorous  development  of  the  abrasion  theory,  a 
copy  of  the  FORTRAN  subroutine  used  to  implement  the  coalescence  formalism,  and 
additional  tabular  and  graphical  data  used  in  the  discussions  in  Chapter  4. 
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Chapter  2 


History  and  Background 


History 

The  first  work  with  high-energy  nucleus-nucleus  collisions  evolved  from  the 
study  of  cosmic  radiation.  Cosmic  radiation  was  conclusively  discovered  in  1912,  when 
Victor  Hess  found  that  ionizing  radiation  levels  actually  increased  with  altitude,  rather 
than  decreasing  as  had  been  expected  if  the  source  was  terrestrial  (Moritz,  1 963).  His 
balloon-based  instruments  determined  that,  at  an  altitude  of  5000  meters,  his  readings 
were  several  times  as  great  as  they  were  at  ground  level.  Initially  thought  to  be  high- 
energy  electromagnetic  radiation,  the  radiation  was  later  found  to  consist  of  highly 
energetic  charged  particles  by  Jacob  Clay.  Early  studies  of  nucleus-nucleus  collisions 
were  made  by  setting  up  targets,  often  in  high  altitude  balloons  or  aircraft,  and  waiting 
for  a  suitable  incident  particle  to  provide  the  reaction  of  interest. 

By  the  1940’s,  early  work  using  cyclotrons  at  Berkeley  and  Yale,  as  well  as 
equipment  associated  with  the  Manhattan  Project  began  to  supplant  the  use  of  cosmic-ray 
targets  with  the  acceleration  of  carbon  and  nitrogen  ions  with  kinetic  energies  as  high  as 
1 0  MeV  per  nucleon.  Later,  more  advanced  accelerators  pushed  the  available  energies 
and  the  mass  of  the  accelerated  ions  to  ever-higher  levels.  Modem  accelerators,  like  the 
Relativistic  Heavy  Ion  Collider  at  Brookhaven  National  Laboratory,  are  capable  of 
accelerating  heavy  nuclei  (such  as  gold)  up  to  the  100  GeV  per  nucleon  range. 
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High-Altitude  and  Space  Radiation 

There  are  three  main  sources  of  radiation  hazards  to  the  crews  of  manned  space  missions 
and  high-altitude  aircraft  flights:  Galactic  cosmic  radiation  (GCR),  solar  particle  events, 
and  charged  particles  trapped  in  the  Van  Allen  Belts.  The  radiations  from  these  various 
sources  share  many  characteristics,  but  each  is  distinct  and  must  be  treated  separately. 
Figure  2.1  shows  the  flux  densities  versus  particle  energy  for  charged  particles  in  the 
space  environment.  Most  of  the  lower  energy  particles  are  protons  or  electrons,  while  at 
very  large  energies  the  particles  can  include  much  heavier  nuclei,  as  discussed  below. 

Ideally,  sufficient  shielding  to  reduce  the  dose  to  reasonable  levels  would  be  used 
to  protect  crews  in  space  and  high-altitude  applications.  Unfortunately,  the  cost  of 
boosting  large  quantities  of  bulk  shielding  into  orbit  or  deep  space  is  far  too  high  to  be 
feasible  (on  the  order  of  $2000/lb).  Hence,  more  accurate  modeling  of  the  radiation 
environment  and  its  interactions  with  spacecraft  structures  will  allow  designers  to 
optimize  the  allocation  of  shielding  materials  and  simultaneously  plan  mission  durations 
that  keep  crew  doses  below  applicable  limits. 

Galactic  Cosmic  Radiation 

Galactic  cosmic  radiation  consists  mainly  of  fully  stripped  ions  (-98%),  with 
energies  up  to  as  high  as  1  trillion  GeV.  The  remaining  portion  (-2%)  is  made  up  of 
electrons  and  positrons.  Of  the  nuclear  portion,  about  87  percent  is  hydrogen,  12  percent 
is  helium,  and  the  remainder  is  heavier  ions.  These  particles  are  believed  to  have  their 
origin  in  astrophysical  phenomena  such  as  novae/supernovae,  neutron  stars  and  black 
holes.  The  relative  abundance  of  these  heavy  ions  drops  fairly  steadily  with  increasing 
mass  up  to  56Fe  (~10‘4x  the  proton  flux),  and  then  drops  off  sharply.  Smaller  spikes  occur 
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-Solar  wind  protons 


-Auroral  electrons 


Trapped  electrons 


100 

Particle  energy,  MeV 


Trapped  protons 
(inner  zone) 


Figure  2.1  The  energy  distribution  of  charged  particles  in  the  space  environment 

(Wilson,  1991) 


with  carbon  and  oxygen  (both  ~10'3x  the  proton  flux).  Solar  activity  and  the  Earth’s 
magnetic  field  serve  to  deflect  the  lower  energy  portion  of  the  flux.  Typical  GCR 
equivalent  doses  in  low  earth  orbit  (for  example,  the  International  Space  Station  or  shuttle 
missions)  are  between  20  and  30  mSv/yr.  For  missions  outside  the  Earth’s  magnetic 
field,  equivalent  doses  may  be  on  the  order  of  40  cSv/yr  at  solar  maximum  and  120 
cSv/yr  at  solar  minimum  (Townsend,  1990). 

Secondary  neutron  production  by  GCR  interactions  with  the  atmosphere  and  with 
aircraft  parts  makes  up  about  half  of  the  total  dose  equivalent  to  flight  crews  and 
passengers  in  high-altitude  aircraft  flights  (NCRP,  1995).  These  doses  have  been  of 
increasing  concern  as  flight  durations  increase  and  the  routes  are  at  higher  altitudes. 
Numerous  studies  have  been  and  are  being  made  to  evaluate  and  predict  the  risk  to  flight 
crews  in  particular  (for  example,  Chee,  2000;  Waters,  2000;  Menzel,  2000). 
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Solar  Particle  Events 


Solar  Particle  Events  (SPE)  are  injections  of  high-energy  particles  into 
interplanetary  space  during  large  disturbances  on  the  sun.  The  frequency  and  intensity  of 
these  events  follow  an  1 1-year  cycle  and  are  rare  during  solar  minimum.  SPEs  are 
usually  associated  with  solar  flares  or  with  coronal  mass  ejections  (CME).  They  are  also 
somewhat  directional  in  nature;  thus,  an  SPE  occurring  on  the  far  side  of  the  sun  is  of 
little  consequence  to  near-Earth  spacecraft  from  a  radiation  protection  perspective.  Most 
SPEs  are  either  too  small  or  contain  too  few  energetic  particles  to  be  of  concern  to  space 
crews.  Extremely  large  SPEs  typically  occur  only  once  or  twice  during  a  given  solar 
cycle.  Near  solar  maximum,  small  SPEs  occur  roughly  hourly,  while  the  largest  occur 
monthly. 

The  spectrum  of  an  SPE  can  vary  quite  a  bit,  but  is  mainly  made  up  of  protons, 
electrons  and  alphas  with  most  particle  energies  below  100  MeV.  Heavier  ions  do  not 
usually  contribute  a  significant  portion  of  the  dose. 

Van  Allen  Belts 

The  Earth  is  partially  shielded  from  external  radiation  by  its  magnetic  field.  One 
side  effect  of  these  fields  is  the  creation  of  several  “belts”  of  trapped  protons  and 
electrons.  Figure  2.2  shows  a  diagram  of  the  particle  distributions.  The  inner  belt  is  the 
more  stable  of  the  two,  and  is  made  up  of  both  high-energy  protons  with  an  average 
energy  of  50  MeV  and  electrons  with  an  average  energy  of  2  MeV.  The  outer  belt  is 
more  diffuse  and  has  mainly  lower  energy  protons  and  electrons.  The  most  energetic 
trapped  particles  are  thought  to  be  the  result  of  the  decay  of  neutrons  produced  by  GCR 
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Distance  from  center  of  Earth,  Earth  radii 


Figure  2.2.  Near-Earth  trapped  radiation  and  solar  proton  environment.  For  clarity,  low- 
energy  integral  fluxes  are  shown  only  on  the  right,  and  high-energy  fluxes  on  the  left. 

(Parker  and  West,  1973) 


interacting  with  the  upper  atmosphere,  while  others  accumulate  from  the  solar  wind,  solar 
particle  events,  and  other  sources. 

Shielding  Considerations 

No  human,  of  course,  travels  unprotected  in  space  or  in  the  upper  atmosphere. 
Vehicle  hulls,  protective  clothing,  equipment,  dedicated  shielding  and  self-shielding  all 
provide  a  measure  of  protection  from  the  radiation  environment.  These  same  shields, 
however,  accentuate  the  need  for  accurate  models  of  neutron  production  in  high-energy 
particle  collisions.  The  typical  shielding  assumed  for  a  spacecraft  without  dedicated 

■y 

radiation  shielding  is  5  g/cm  of  aluminum.  Plans  for  long-duration  space  missions 

<y 

include  provision  of  a  “storm  shelter”  with  1 8  to  24  g/cm  of  aluminum  for  use  during 
SPEs  (Simonsen,  1991). 
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Chapter  3 

Neutron  Production  Model 


Particle  production  in  high-energy  nucleus-nucleus  collisions  is  considered  to 
come  from  two  stages:  first,  nucleons  are  sheared  off  (abraded)  during  collision,  then  the 
remaining,  excited  pieces  of  the  projectile  and  target  decay  by  particle  and  gamma 
emission. 

The  present  model  includes  the  abrasion  (knockout)  step,  the  ablation  step,  and 
then  adds  a  coalescence  step.  Previous  versions  of  these  models  (Cucinotta,  1997) 
applied  the  abrasion  (knockout)  portion  to  both  the  projectile  and  the  target  and  ablation 
to  the  projectile  only,  but  did  not  incorporate  a  coalescence  step.  Incorporating 
coalescence  is  the  major  new  feature  of  this  work.  In  the  abrasion  stage,  a  Glauber  model 
is  used  that  incorporates  effects  of  final-state  interaction  (Glauber,  1970).  A  more 
detailed  derivation  of  the  abrasion  model  due  to  Cucinotta  (Cucinotta,  1 997)  appears  as 
Appendix  A.  The  ablation  stage  uses  a  classical  evaporation  model  (Cucinotta,  1994a; 
Weisskopf,  1 940)  to  evaluate  the  decay  of  the  remaining  portions  (called  prefragments) 
of  the  target  and  projectile  nuclei.  Finally,  the  coalescence  stage  uses  a  modified  version 
of  an  existing  model  (Nagamiya,  1981)  to  include  particles  emitted  by  the  projectile  as 
well  as  the  target. 

Abrasion  Theory 

In  an  abrasion-ablation  model  the  momentum  distribution  of  the  nucleons 
produced  by  the  collision  is  given  by 
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In  the  Glauber  model  the  scattering  amplitude  for  knockouts  is  related  to  the 
phase  space  of  each  particle  appearing  in  the  final  state.  Neutrons  produced  by  this 
model  are  considered  to  originate  in  the  projectile  rest  frame;  by  exchanging  the  target 
and  projectile  and  converting  to  the  target  rest  frame  (lab  frame),  values  for  the  target  can 
be  obtained.  Using  a  single  momentum  vector  denoted  as  px  to  apply  closure  to  the  target 
final  state,  the  cross  section  is  given  by 


f  fe  -  e,  >  (p,  -  P/  k/  (2) 

P  X  »= 1  j=\ 

where  P  is  the  relative  projectile  velocity  divided  by  the  speed  of  light,  F*  represents  the 
prefragments,  n  is  the  number  of  nucleons  knocked  out  of  the  projectile  in  the  collision 
with  the  target,  T  is  the  transition  amplitude  that  describes  the  collision,  and  i  and  f 
denote  the  initial  and  final  states,  respectively.  Using  a  Glauber  model  for  the  transition 
amplitude  (T)  in  equation  (2),  along  with  closure  on  the  target  states,  allows  the 
momentum  distribution  from  abrasion  to  be  written  as 


where  Ap  is  the  projectile  mass  number,  PA,'~n(b)  denotes  the  probability  that  Ap  -  n 
nucleons  are  not  removed  by  knockout  at  impact  parameter  b,  and  the  distribution  of 
knockouts  versus  impact  parameter  is  given  by  dNn/dk.  The  latter  can  be  written  in  terms 
of  single  collision  terms  dNi/dk  as 


dK 

dk 


dA\ 

dk 


[i-p(b)r 


(4) 
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In  the  optical  limit,  we  have 


P(b)  =  exp[-2Im*(b)l  (5) 

L  ap 

where  the  eikonal  phase  x  is  given  in  terms  of  one-body  form  factors  F(q)  for  the 
projectile  and  target  and  the  nucleon  scattering  amplitude  fjsrNi  as 

x<P)  =  7—— ;  J^cxpO'q-h  WpiaWMfmtiq)  (6) 

(2 nkNN)  J 

with  At  as  the  target  mass  number  and  q  denoting  the  momentum  transfer  in  the  collision. 
The  form  factor  F(q)  is  the  Fourier  transform  of  the  appropriate  (projectile  or  target) 
nuclear  density. 

Ablation  Theory 

After  the  projectile-target  collision,  a  distribution  of  compound  (prefragment) 
nuclei  in  various  states  of  excitation  remains.  The  excitation  spectrum  is  treated  as  an 
average  state  in  the  abrasion-ablation  model,  with  a  single  or  small  number  of  average 
excitation  energies  used  to  describe  the  prefragment  and  the  strength  of  the  state 
determined  by  the  total  abrasion  cross  section.  Next,  estimate  the  average  excitation 
energies  and  calculate  the  neutron  spectrum  from  ablation  with  the  Weisskopf-Ewing 
statistical  decay  model.  Define  Pn(j,  k)  as  the  probability  that  a  prefragment  labeled  j 
with  mass  number  Aj,  charge  number  Zj,  and  excitation  energy  Ej  ,  emits  a  neutron  of 
momentum  k  (Weisskopf,  1940;  Kikuchi,  1968).  The  momentum  distribution  for  neutron 
production  from  ablation  in  the  projectile  rest  frame  is  written  as: 

-7-  =’Z^,(AJ,ZJ,E])P,Xj,k)  (7) 

V  «k  )  abl  j 
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The  total  abrasion  cross  section  in  equation  (7)  is  evaluated  from  the  abrasion  momentum 
distribution  in  equation  (3). 

In  the  statistical  model,  the  prefragment  (compound  nucleus)  is  assumed  to  be 
infinitely  heavy,  and  the  emission  spectrum  is  assumed  to  be  isotropic  in  the  rest  frame. 
The  probability  function,  as  given  in  (Weisskopf,  1940;  Kikuchi,  1968),  is 


Pn(j\E„)  = 


2M„g„E,^cNwo(E]  -E„) 


(8) 


where  pn  is  the  reduced  neutron  mass,  gn  is  the  statistical  weight,  ctcn  is  the  formation 
cross  section  by  the  inverse  process,  wQ  is  the  level  density  of  the  residual  nucleus,  and 


F,  =  f'S‘P,U,E)dE 


(9) 


Competition  occurs  between  the  emission  of  the  light  particles  n,p,  2H,  3H,  3 He, 
and  a.  The  model  parameters  are  taken  from  (Dostrovsky,  1959).  If  sufficient  excitation 
energy  is  available,  then  several  neutrons  will  be  emitted.  The  cumulative  spectrum  is 


(10) 


In  this  work,  terms  through  the  third  order  in  the  series  of  equation  (10)  are  considered. 
The  neutron  momentum  distribution  in  the  lab  frame  is  found  by  multiplying  equation  (7) 
by  the  appropriate  neutron  energy  to  form  a  Lorentz  invariant  and  then  performing  the 
transformation  to  the  laboratory  system. 

Coalescence  Theory 

The  abrasion  and  ablation  models  presented  thus  far  yield  single  particles  with  a 
forward-peaked  distribution.  Certain  of  these  particles,  with  similar  momenta,  may 
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coalesce  to  form  heavier  nuclei,  such  as  deuterons,  tritons,  helions,  and  alphas. 
Schwartzschild  and  Zupancic  (Schwatzschild,  1963)  predicted  that  the  deuteron  density 
in  momentum  space  is  proportional  to  the  proton  density  times  the  probability  of  finding 
a  neutron  within  a  small  sphere  of  radius  po  around  the  proton  momentum.  Following 
Awes  (Awes,  1981),  the  probability  P  of  finding  a  single  nucleon  in  the  coalescence 
volume  centered  at  a  momentum  per  nucleon  p  is  given  by  the  product  of  this  volume 

4tf  3 

( —  p0  )  with  the  single  nucleon  momentum  density 


_  4tt  3  Y  d3N(p) 

3  Po  m  dp 3 


(11) 


where  d3N(p)/dp3  represents  the  differential  nucleon  multiplicity,  m  is  the  average 
nucleon  multiplicity,  and  y  =  1  +  77  m.  For  y  only,  m  is  the  nucleon  mass,  rather  than  the 
multiplicity.  This  is  simply  the  volume  times  the  number  of  nucleons  per  unit  volume,  all 
in  momentum  space. 

For  a  given  multiplicity  m,  such  as  when  m  nucleons  are  produced  in  an  event,  the 
probability  of  finding  n  of  them  ( n  <  m)  in  the  coalescence  volume  will  be  given  by  the 
binomial  distribution 


P(n  |  m) 


f m ^ 


p"(\-py 


(12) 


\n  J 

Each  multiplicity  will  have  a  probability  f(m)  of  occurrence,  which  can  be  summed  over 
the  distribution  of  multiplicities  to  obtain  the  mean  probability  for  finding  n  nucleons  in 
the  coalescence  volume. 
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In  the  case  of  low  average  multiplicities,  it  is  reasonable  to  assume  a  Poisson 
distribution  of  multiplicities 

,/  \  -s 

m\ 

which  can  be  substituted  into  equation  (13)  which  gives 

-  V  (w)  "  m- 


(14) 


m\  n\(m-n). 


Pn(\-Pf~" 


n\  7>qV\  ’  J 


(15) 


_  [mP]  ( 


-mP 


n\ 

which  is  exact  for  a  Poisson  multiplicity  distribution.  At  relativistic  energies  large 
multiplicities  are  observed  with  a  non-Poisson  distribution,  in  which  case  it  is  assumed 
that  equation  (13)  can  be  approximated  by 


./  \ 

m 

\n\  m)=\ 

P"(l 


(16) 


Then  if  P  « 1  and  m»n  the  binomial  distribution  becomes  Poisson  in  form,  which 
again  results  in  equation  (15). 

Since  mP  is  quite  small,  the  exponential  term  in  equation  (15)  can  be  ignored, 


which  means  that  the  average  probability  for  having  N  neutrons  and  Z  protons  in  the 
coalescence  sphere  is 


(17) 


(/.(Ar,Z))=kiIk5J' 


Z!  N\ 

where  it  is  assumed  that  the  probabilities  for  the  presence  of  neutrons  and  protons  are 
independent.  In  this  context,  P(N,  Z)  represents  the  probability  of  forming  a  composite 
particle  with  the  momentum  per  nucleon  p. 

Substituting  equation  (11)  into  equation  (17)  and  dividing  by  the  coalescence 
volume  gives  the  composite  particle  momentum  distribution  in  the  form  of  the 
coalescence  relation  used  at  relativistic  energies 


d2N(Z,N ) 
dp 3 


fZ.+zAZ  1  (An 


\N'+NPj 


N\Z\ 


PoY 


\A~V  ^ 


d3N(0,l) 
dp 3 


(18) 


where  Z,  N,  and  A  are  the  composite  particle  charge,  neutron  number  and  mass 
respectively;  Np  and  N,  are  the  projectile  and  target  neutron  numbers,  and  Zp  and  Z,  are 
the  projectile  and  target  charge  numbers.  The  multiplicity  is  related  to  the  momentum 
distributions  by 

d3N  _  1  d3cr 
dp 3  cr0  dp 3 

The  quantity  o0  is  the  projectile-target  total  reaction  cross  section,  and  is  calculated  using 
one  of  two  formulations,  described  below.  Applying  equation  (19),  we  obtain 


(20) 

Since,  as  part  of  the  abrasion  and  ablation  steps,  we  obtain  Lorentz  invariant  cross 
sections  for  the  total  neutron  production,  we  want  to  convert  this  to  utilize  a,  (n) ,  the 


£/3<t(Z,A) 

f  z,+zp} 

z 

1 

(a  3  > 

4^o  Y 

A- 1 

f  ^3cr(0,l)^ 

dp 3 

{ N’+NP) 

N\Z\ 

K  3CJ0  J 

l  ¥  J 
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invariant  neutron  cross  section.  The  invariant  cross  section  is  related  to  the  triple 


differential  cross  section  by  =  —  <r, .  Substituting  into  equation  (20),  we  get 

dp  E 


4r^(A)  = 


z,+zp 

kN‘+Np; 


N\Z\ 


/  3 

1  4  7?p0  y  '  ' 

V  3cto  j 


\A 


(21) 


(J,(A)  = 


(22) 


but  Ea  =  AE : 


EA  _EA  _  A 


Ea  Ea  Ea~] 


so 


(7,  (A) 


ACT 

ea~\ 


(23) 


which  yields 


cr,{Z,N) 


( Z'+V 

z 

A 

"4^?0V 

A-\ 

N\Z\ 

v  3ct0£  , 

(v,(nj)A 


_ m_ 

T  +  m 


— ,  the  inverse  of  the  nucleon  rest  mass  energy.  Thus, 
m 


(24) 


a,(Z,N) 


z,+zp 

z 

A 

(  4^03 1 

u+ivJ 

N\Z\ 

3(7 

(<7,{n))A 


(25) 


Nagamiya  et.  al.  rewrote  equation  (25)  in  a  slightly  different  form  by  using  a 
constant,  energy-independent  coalescence  coefficient  Ca  in  an  equation  of  the  form 


(  ^ 
a  cxA 

=cA 

f 

E 

\  dP  A  J 

K. 

)) 

(26) 
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where  pA  =  App  and  A  is  the  mass  of  the  coalesced  fragment.  Comparing  equation  (26) 
with  equation  (21),  we  can  see  that  CA  is  comparable  to  our  CT,  which  is  also  energy  and 
angle  independent. 

Parameterization  of  oo 

For  nucleon-nucleus  reactions,  a  parameterization  by  Wilson  (Wilson,  1988)  was 

used, 

o-0  =45/(£)^°7{l  +  0.016sin[5.3-2.631n(^)]}  (27) 

where  A  is  the  mass  number  of  the  target  nucleus  and  the  energy-dependent  factor  is 

/(£)  =  1-0.62  exp(-£  /  200)  sin(l  0.9£'° 28 )  (28) 

with  the  incident  kinetic  energy  E  in  units  of  MeV/nucleon. 

For  nucleus-nucleus  reactions,  a  different  parameterization  was  used  (Townsend, 

1986), 

a0  =  nr3  f3{E)[A)i3  +  A]'3  -  S( A p ,  Ar ,  E)f  (29) 

with 


J3(E)  =  l  +  5E- 


(30) 


and 


6  i 


S(Ap,AT,E)  =  0.200  + 


\apJ 


f  i  ^ 


+ 


\At  j 


-  0.292  exp  (-E  /  792)  cos(0.229 E04  53 ) 


(31) 


with  r0  =  1.26  fm  and  E  again  in  units  of  MeV/nucleon. 
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Chapter  4 

Results  and  Discussion 


In  order  to  determine  the  accuracy  and  validity  of  this  neutron  production  model, 
a  comparison  with  experimental  data  is  conducted  for  four  interactions:  357  MeV/A 
Ca+H  (Tuve,  1997),  565  MeV/A  Ca+H  (Tuve,  1997),  75  MeV/A  C+Al  (Heilbronn, 
unpublished),  and  390  MeV/A  Ne+Ne  (Madey,  1985). 

Model  Implementation  Considerations 

The  implementation  of  the  complete  model  consists  of  about  2800  lines  of 
FORTRAN  77  code,  which  was  developed  and  run  on  a  DEC  Personal  Workstation  433a, 
compiled  with  Compaq  Fortran  for  Alpha  Linux  vl  .0  and  run  on  RedHat  Linux  v7.0. 
Results  were  transferred  to  Microsoft  Excel  2000  for  plotting  and  statistical  analysis.  The 
subroutine  calculating  coalescence,  the  major  contribution  to  this  model  in  this  work,  is 
included  in  appendix  B. 

As  with  any  computer  model  of  such  complexity,  some  numerical  problems  can 
creep  in  to  the  code.  Roundoff  error  is  the  most  common  trouble.  For  example,  at  one 
place  in  the  ablation  code,  the  numerator  and  denominator  of  a  particular  expression  were 
supposed  to  proceed  “in  step”  to  0.  L’Hopital’s  rule  analytically  showed  the  limit  to  be 
0;  however,  roundoff  error  generated  a  result  where  some  of  the  cross  sections  at  small 
angles  behaved  like  lim(l  /  x) .  This  was  corrected  by  numerical  interpolation  from  the 

jt-»0 

surrounding  cross  section  estimates. 


17 


Another  numerical  problem  was  that  the  code  that  modeled  the  abrasion  behavior 
of  the  target  produced  results  which  behaved  asymptotically  at  large  lab  angles  and 
eventually  dominated  the  overall  spectrum.  This  problem  has  not  been  corrected  at 
present  and  remains  for  future  work. 

As  discussed  in  the  previous  chapter,  the  coalescence  step  of  the  present  model 
depends  on  the  selection  of  a  single  parameter,  po,  the  coalescence  radius.  A  wide  range 
of  values  for  po  has  been  proposed;  some  authors  propose  different  values  for  each 
composite  particle  formed.  Awes  et  al.,  for  example,  proposed  values  of  170  MeV/c,  215 
MeV/c,  and  270  MeV/c  for  deuterons,  tritons,  and  alphas  respectively  (Awes,  1981). 
Gutbrod  et  al.  proposed  values  ranging  from  106  to  147  MeV/c  for  a  range  of  interactions 
(Gutbrod,  1976).  Nagamiya  et  al.,  whose  results  were  most  influential  on  the  present 
work,  selected  a  single  value  of  90  MeV/c  and  obtained  very  good  agreement  with  their 
data  (Nagamiya,  1981).  For  the  purposes  of  comparison,  results  here  are  presented  with 
values  for  po  of  90  MeV/c  and  125  MeV/c,  the  latter  being  representative  of  the  values 
found  by  Gutbrod  et  al. 

Results  for  Calcium-Hydrogen 

Figures  4.1  (a  and  b)  and  4.2  (a,  b,  and  c)  display  the  results  obtained  for  357  and 
565  MeV  per  nucleon  calcium  projectiles  on  hydrogen  targets  integrated  out  to  3°  from 
the  beam  angle.  The  figures  also  include  the  results  obtained  with  the  earlier  version 
(Cucinotta,  1997)  of  the  current  model,  prior  to  the  inclusion  of  coalescence.  As  these 
figures  demonstrate,  the  current  work  fits  the  high  energy  portion  of  the  spectrum  much 
more  closely  than  the  previous  model,  although  the  high-energy  portion  still  tends  to  be 
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c)  Angular  distribution  results 
Figure  4.2  continued 
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Table  4.1  Chi-squared  goodness  of  fit  check  for  previous  and  current  models.  Bold  type 
indicates  best  fit. 


Beam  Energy 

t 

xz 

T  (MeV  per  nucleon) 

E<T 

E>T 

All  E 

357  Previous  Model 

4.91 

128.79 

133.70 

357  This  Work,  p0=  125  MeV/c 

17.33 

5.90 

23.23 

357  This  Work,  p0  =  90  MeV/c 

9.86 

8.39 

18.25 

565  Previous  Model 

769.0 

8449.0 

9218.0 

565  This  Work,  p0=  125  MeV/c 

1759.4 

340.2 

2099.6 

565  This  Work ,p0  =  90  MeV/c 

1491.9 

491.0 

1982.9 
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overestimated.  Since  prior  model  data  was  available  for  these  projectile-target 
combinations,  a  chi-squared  test  was  made,  with  the  results  in  table  4.1  below. 

Both  the  figures  and  the  table  above  demonstrate  that,  for  both  of  these  projectile- 
target  combinations,  using  a po  of  90  MeV/c  is  a  better  choice  than  125  MeV/c.  Both 
choices  were  a  better  fit  to  the  energy  spectrum  above  the  beam  energy  than  was  the 
previous  model,  but  the  previous  model  was  a  better  fit  below  the  beam  energy.  Figure 
4.2c,  the  angular  distribution,  shows  that  none  of  the  models  predicts  a  sufficient  falloff 
with  increasing  angle.  It  should  be  noted  that,  for  the  357  MeV  per  nucleon,  there  were 
only  two  experimental  data  points  below  the  beam  energy,  making  the  comparison 
somewhat  weak.  Overall,  the  current  model  with  a  po  of  90  MeV/c  is  the  best  fit  at  both 
357  MeV  per  nucleon  and  at  565  MeV  per  nucleon.  Better  low-energy  data  for  these 
projectile-target  combinations  would  help  clarify  the  situation. 

Results  for  Carbon-Aluminum 

Figures  4.3,  4.4  and  4.5  compare  the  results  for  75  MeV  per  nucleon  carbon 
projectiles  on  aluminum  targets  at  10,  60  and  125  degrees.  Appendix  C  contains 
additional  figures  for  30,  45,  and  90  degrees.  Several  aspects  of  these  figures  warrant 
discussion.  As  the  experimental  data  are  recent  and  unpublished  (Heilbronn,  private 
communication),  they  have  been  arbitrarily  scaled  so  that,  at  each  lab  angle,  the 
maximum  experimental  value  is  1.0.  Also,  because  the  experimental  data  analyses  are 
not  yet  complete,  the  final  values  may  change.  Figures  4.4,  and,  in  particular,  4.5 
demonstrate  the  problem  of  asymptotic  behavior  in  the  target  knockout  portion  of  the 
code. 

A  statistical  analysis  comparing  the  y2  values  for  po  =  90  MeV/c  and  for 
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Figure  4.3  Comparison  of  results  for  75  MeV  per  nucelon  C+Al  at  10“  lab  angle 


■  Experiment 

+  This  work  (Po=125  MeV/c) 
-  This  work  (Po=90  MeV/c) 


1.00E+00 


+  +  +  +  -*-  +  + 


’  ’  n?r~— ,, 

i  ,  — - 


20  40  60  80  100  120  140  160 

E,  MeV 


Figure  4.4  Comparison  of  results  for  75  MeV  per  nucleon  C+Al  at  60  lab  angle 


Figure  4.5  Comparison  of  results  for  75  MeV  per  nucleon  C+Al  at  125°  lab  angle, 
showing  the  asymptotic  behavior  of  the  present  model  at  large  lab  angles. 


p0=  125  MeV/c  (included  in  appendix  C)  shows  that  125  gives  the  better  fit  for  all  but 
the  125°  data;  however,  above  60  degrees  the  fit  is  too  poor  to  be  significant. 

Results  for  Neon-Sodium  Fluoride 

Figure  4.6  shows  the  results  obtained  for  390  MeV  per  nucleon  neon  projectiles 
on  sodium  fluoride  targets  at  a  0°  lab  angle.  For  this  model,  experimental  data  for  Ne- 
NaF  was  compared  with  computed  results  for  Ne-Ne,  which  is  nearly  the  same  for 
analysis  purposes.  Because  numerical  values  for  the  experimental  data  uncertainties  were 
not  available,  no  formal  statistical  analysis  was  made.  However,  it  is  clear  from  the  plot 
that  the  selection  of po  =  90  MeV/c  is  the  better  choice  for  this  projectile-target 
combination.  In  this  case,  as  well,  the  underestimation  at  the  low  energies  is  quite 
evident,  as  is  the  smaller  overestimation  at  high  energies. 
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Chapter  5 


Conclusions  and  Future  Work 


Conclusions 

The  expanded  model  described  in  the  preceding  chapters  provides  a  fairly  good  fit 
to  published  data,  although  further  development  is  needed.  In  most  of  the  cases 
considered,  the  addition  of  coalescence  and  the  target  ablation  provided  a  better  fit  than 
the  previous  model,  with  the  selection  of  90  MeV/c  for  the  coalescence  radius  generally 
yielding  the  best  results.  The  good  agreement  between  the  model  and  published 
experimental  data  suggests  that  significant  progress  has  been  made  in  developing  an 
accurate  tool  to  predict  neutron  spectra  from  heavy  ion  collisions. 

Future  Work 

Three  code-specific  problems  need  to  be  addressed.  The  first,  already  mentioned 
in  chapter  4,  is  to  provide  a  better  solution  to  the  problem  of  the  cross  section  becoming 
infinite  at  energies  very  close  to  the  beam  energy  for  very  small  lab  angles.  The 
linearization  already  implemented  in  the  code  should  either  be  replaced  by  an  appropriate 
curve  fit  or  a  better  approximation  for  the  analytical  expression  should  be  found. 

Second,  a  better  formulation  for  the  target  abrasion  model  should  be  developed. 
The  present  code  yields  cross  sections  that  behave  asymptotically  at  large  lab  angles, 
whereas  the  cross  sections  should  be  dropping  to  0  at  higher  energies.  The  values 
involved  in  the  calculations  are  small  enough  that  roundoff  error  may  be  involved. 

Third,  due  to  inherent  limitations  in  the  way  FORTRAN  77  assigns  array  sizes, 
the  code  is  currently  restricted  to  a  maximum  nucleus  mass  (target  or  projectile)  of  68. 
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Recasting  the  code  into  a  newer  version  of  FORTRAN  would  enable  array  sizes  to  be 
dynamically  assigned,  and  would  remove  this  limitation. 

Aside  from  code  modifications,  further  research  is  needed  to  select  a  better  value 
for  the  coalescence  radius  po,  or  a  parameterization  if  a  single  value  is  not  appropriate. 
As  one  might  expect,  the  larger  the  value  chosen  the  larger  the  number  of  nucleons  that 
coalesce  to  product  light  nuclei.  For  most  of  the  cases  considered  herein,  p0  =  90  MeV/c 
gave  a  good  fit  to  the  experimental  data,  but  not  in  every  case.  The  values  in  the 
literature  vary  greatly  with  projectile  and  target  masses,  beam  energies,  and  type  of 
composite  fragment  formed. 

Also  desirable  would  be  more  detailed  data  concerning  the  low  energy  neutron 
spectrum  to  help  determine  the  accuracy  of  the  model  in  this  region. 

Once  a  fully  acceptable  model  has  been  produced,  a  database  of  neutron 
production  cross  sections  for  a  wide  range  of  projectile-target  and  energy  combinations 
can  be  developed  for  use  in  accurately  predicting  doses  from  cosmic  radiation  and  large 
solar  particle  events.  The  model  can  also  be  used  to  predict  light  ion  production  cross 
sections  produced  by  the  coalescence  step. 
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Appendices 


Appendix  A 


Abrasion  Theory 


(From  Cucinotta,  1997) 


In  the  Glauber  model  the  scattering  operator  for  nucleus-nucleus  collisions  is 
(Glauber  1970,  Cucinotta  1994b) 

/=—  frf2*exp(/q.b)r(b)  (Al) 

where  K  is  the  relative  wave  number  of  the  projectile  target,  b  is  the  impact  parameter, 
and  q  is  the  momentum  transfer.  The  profile  function  representing  the  multiple¬ 
scattering  series  at  high  energies  is 

r(b) = i  -  P[P  -  rfJ  (■>  -  «„  - »,)]  (A2) 

C,i 

where  £  and  j  label  the  target  and  projectile  constituents,  respectively.  In  equation  (A2), 
TSJ  is  the  two-body  profile  function  with  the  internal  coordinate  having  components 
r=(s,  z)  and  represents  individual  nucleon-nucleon  collisions  occurring  during  the 
interaction  of  the  nuclei. 

The  scattering  amplitude  of  equation  (Al)  is  related  to  the  production  cross 
section  for  a  projectile  nucleon  from  the  abrasion  (knockout)  process  by 

—  =  y  — \dE.d2qd2bd2b'exp[iq  •  (b  -  b’)]J(£,  -  E, ) 
dk  Y(2 nf  J  F  '  } 

\(rp\r'  (b’  )\XF'kj  ^(kjF'X |T(b)|Pr)  (A3) 

where  kj  denotes  the  wave  numbers  of  the  abraded  nucleons,  F*  denotes  the  prefragment 
(with  Af*  =  Ap  -  n),  and  initial  and  final  states  have  been  inserted  into  equation  (A3). 
Next,  the  change  in  energy  of  the  target  from  the  collision  is  as  follows: 
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Er-E„  =  Et  -V(pT-q )J+M}  =  ET  1  ~y+Ml  M * 

V  hr 


where  Mj  and  Mx  are  the  target  masses  for  the  initial  and  final  state,  respectively.  After 
closure  on  the  final  target  states, 

~  =  \dEFJ2qd2bd2b'Qxp[i(\\b-b')}(Tn^,b' ,k,q,EFt)  (A5 

dk  (2n)  J 

where 


"  f  dk: 


cr„(b,b',k,q,£/,)  =  (r||  fn  ^  -^)(^|rUb')|F\y.)(k,F’|r(b)|p)[|r) 


To  consider  the  energy  conservation  term  8 (ErEj)  in  equation  (6)  above,  use  the 


Fourier  transform  pair 


cr„(0=  §dEexp(-iEt)an(E) 


v„(E)  =  J ~  exp(fEY  )<x„  (t ) 

J  ?  rr 


In  the  projectile  rest  frame,  we  have 


«  k^  _  k^ 
E,-E/=S.-Tf.+Er-Er-X^  =  E-X^ 

7=1  j=wz-mN 


where  S„  is  the  separation  energy  and  7>*  is  the  recoil  energy  of  the  prefragment 
(including  excitation  energy). 


From  Eqs.  (A6),  (A7)  and  (A9), 
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-,.(')=<i  jn 


7=2  l 


dkj 

(2^ 


exp(-fZ 


k]t 


M2  mN 


r+ (b1  )\f\j  YjljF*  |r  (b)|p)||r )  (A  1  o) 


To  simplify,  first  factor  the  profile  function  into  projectile  participant  and 


spectator  terms  as 


r(b)=i-n&(b-s;)fle, (•>-»,) 


(All) 


l=n+\ 


7=1 


where 


ft=f[0-r„)  (M2) 

<r=i 

where  Q(b  -  s)  is  optically  averaged  over  target  coordinates. 

If  the  prefragment  mass  is  much  greater  than  that  of  the  knockouts,  the  projectile 
wave  function  can  be  factored  as 


(A  13) 


where  |  F)  and  |  <j>n )  are  the  wave  functions  of  the  spectators  (core)  and  participants 

(knockouts)  respectively.  Using  plane  wave  states  for  |  k ;  ^  and  substituting  equations 
(A1 1)  and  (A13)  into  equation  (A10),  we  obtain 


efrjrfrj ’exp(/k  ■  Xj )qJ (b'-s^X^b -Sj) 


x  n 

7=2 


elk 


(2tt) 


3  dr j  dr j'exp(ik  j  xy  )exp| 


it/-,,, 


<t, 


P'j  (b'-s  J ' )Qj  (b'-s  j ')  yn  (r, rn ' )<j>n  (rj , . . . ,  r„ ) 

(A  14) 
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where  \j  =  rj  -  rf.  Using  the  coherent  approximation  for  the  target  wave  function  in  the 
intermediate  states  and  the  independent  particle  model  for  the  fragment  wave  function 
leads  to  the  following  equation: 

'A  A  dN 

an(t)=  '  *P A’’~" (b, b' ) Ah_, (b, b' , /)  — 1  (A15) 

ynj  dk 

where  the  function  Piby)  describes  the  projectile  spectators  as  given  by 

|nao.-s,)|^>  <ai6> 

/  i 

After  performing  closure  in  equation  (A16),  we  find  that 

!P^(b,b-)  =  (rF|n^(b’-s/)na(b-s/)|^>  (A17) 

i  i 

In  equation  (A15)  we  have  defined 

=  f**'exp(ik  x)p(r,r')e,,(b'-s')a(b-s)  (A  18) 

dk  {2n)  J 

where  p(r  ,r’)  =  <p*(r’)0(r)  is  the  one-body  density  matrix  of  the  projectile.  Next,  from 
equation  (A  14),  after  evaluation  of  the  integrals  over  k7  for  j  >2,  and  transforming  back  to 
energy  space,  we  find  that 

A„_,  (b,  b' ,  7>„  k)  =  (T  |  JfJ  dtjdtj '  p(r;  ,r/)  Q)  (b'-s } ' ) 

j= 2 
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and 


Vi  =,  Ex7 

>2 


(A21) 


and  where  is  the  cylindrical  Bessel  function  of  the  first  kind  of  order  m. 

For  n=l,  we  have  A0  =  S  ( 7> *).  If  we  assume  forward-peaked  density  matrices  (about),  a 
small  argument  expansion  of  the  Bessel  functions  can  be  developed  (#ref)  that  results  in 


A,,-!  (b,  b' ,  Tj?t ,  k) «  C„_ 


Trt+- - Sn-(Er-Ex) 

2  mN 


f 


K 


n- 1 


4-, 


\ 


b,b', 

v  v«- 1 


+0(£V) 


(A22) 


7T  ^3"  TC 

where,  for  example,  C\  =  1 ,  Cj  =  — ,  C,  = - ,  and  C4  = - . 
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The  nucleon  momentum  distribution  from  abrasion  then  takes  the  following  form: 


v  dk  Jabr  „ 


v  n  j 


{In') 


dN | 

dk 


(A23) 


If 

q1  +  M\  -  M\  «E). 

then  equation  (A23)  can  be  approximated  by 
^  dcr^ 


(A24) 


V  dk  J  abr  "  V 

To  include  the  final-state  interaction  of  the  nucleon  knockouts,  use  the  Eikonal 
model  described  elsewhere  (Cucinotta  1994b).  In  this  model,  the  plane  waves  are 
replaced  by  the  distorted  wave  for  the  nucleon-projectile  recoil  interaction  evaluated  at 


=  Z  t  |K^“"(b)^  fAA-A^k) 


(A25) 
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the  relative  energy  between  the  knockout  and  the  recoil.  Modifying  equation  (A1 8)  as  in 
reference  (Cucinotta,  1 994b)  gives 

^  =  [drdr'exp(ik-x)p(r,r')exp[-2lmx(-)(y)]Ql  (b'-s')Q,(b-s)  (A26) 

ak  (2ft)  J 

where  is  the  outgoing  Eikonal  phase.  Equation  (A26)  ignores  off-shell  effects  but 
includes  the  energy  dependence  of  the  final-state  interaction  and  assumes  a  medium 
modified  interaction,  as  described  in  reference  (Cucinotta,  1994b). 

The  optical  limit  of  the  profile  functions  may  be  used  for  ApAj»\  (Glauber, 

1970;  Cucinotta,  1992),  which  gives  us 

P+  (b,  b')  =  exp{i[Z(b)  -  /  (b1)]}  (A27) 

where  the  Eikonal  phase  is 


where  crNN  is  the  two-body  total  cross  section,  km  is  the  relative  momentum  in  the  two- 
body  center  of  mass  frame,  B  is  the  two-body  slope  parameter,  and  a  the  ratio  of  the  real 
to  imaginary  parts  of /v,v(q  =  0).  For  inelastic  terms,  we  write  in  the  optical  limit 
(Cucinotta,  1992)  that 

Q]  (b'-s/)S,  (b  -  »J )  =  exp[Q(b'-s/ ,  b  -  s ,  )|  ■ - 1  (A30) 

with 
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a(b'-s  ',b-s  ■) 


1  2  £  f  d2qd2q'exp[iq  ■  (b  -  Sj  +  sa )]  exp[-/q  •  (b’-s  j  '+sa )]  fJjN  (q'  )fNN  (q) 


(2 *kNNr  « 


(A31) 

The  one-body  density  matrix  is  evaluated  using  a  local  density  approximation, 
which  is  assumed  to  factor  as: 

P(r,r>p(y)n(x)  (A32) 

with  x  =  r  -  r’  and  y  =  ,5(r  +  r’)  and  where  the  one-body  density  is  given  by  the  diagonal 
part  of  the  density  matrix 

P(r)  ~  p(r,r'),  (r'=r)  (A33) 

and  n(\)  is  the  Fourier  transform  of  the  nucleon  momentum  distribution 

n(x)  =  [(ip  exp(/p  ■  x)n(p)  (A3  4) 

where  «(p)  is  defined  by 

«(p)  =  Jjrc/r'  exp(-/p  •  x)  p(r,  r' )  (A3  5) 

with  normalization 


(A36) 

The  one-body  density  and  the  nucleon  momentum  distribution  at  small  to  modest  values 
of  p  are  known  from  electron  scattering.  For  large  values  of p,  backward  production  of 
protons  suggests  that  large  enhancements  to  the  nucleon  momentum  distribution  result 
from  correlation  effects  when  compared  with  predictions  of  independent-particle  models. 
In  addition,  this  enhancement  is  largely  independent  of  nuclear  mass. 

Haneishi  and  Fujita  (Haneishi,  1986)  have  introduced  the  following  momentum 
distribution: 
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(A37) 


3 

«(P)  =  «oZC'exP 

#=1 

where  no  is  a  normalization  constant.  The  last  term  in  equation  (A3  7)  is  expected  to 
directly  reflect  the  nuclear  correlations.  The  value  of  C3  =  0.003  used  in  (Haneishi,  1986) 
was  found  to  be  too  large  in  (Cucinotta,  1994a;  Cucinotta,  1994b);  thus  the  value  of  C3  in 
(Cucinotta,  1997)  was  set  at  0.0008. 
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Appendix  B 

Coalescence  Subroutine 


The  following  FORTRAN  77  subroutine  was  used  in  the  determination  of  the 
coalescence-reduced  neutron  production  cross  sections  discussed  in  chapter  4.  Please 
note  that  the  real*  16  type  declaration  and  the  qsqrt  function  call  may  not  be  available 
all  FORTRAN  77  implementations. 


c*  *  SUBROUTINE  COAL(iat,iap,izt,izp,adeg,tnnlab,ds_total,linvt,eel,ne)*  * 
c  this  subroutine  calculates  a  reduction  of  the  neutron  production  cross 
c  sections  due  to  coalescence  based  on  equations  by  Awes,  Nagamiya,  and  Gutbrod 
c  returns:  ds_total 
c  called  by:  ubernspec 
c  calls:  none 

SUBROUTINE  COAL(iat,iap,izt,izp,adeg,tnnlab,ds_total,linvt,eel,ne) 
implicit  none 

real*8  pi, r0,rf2,sigma0,eel(50),f, beta, delta, nucleonmass, 

&proton_factor,p0,adeg,plab(50),temp_total,dstotal(50),linvt(50), 

&gamma,tnnlab 

real*  1 6  k(4),rfl ,r,s,t,u,q,n_obsl ,n_obs2,eterm,elab,p_obsl ,terml , 
&term2,term3, root, parta,p_total,p_total2,partb, inner 
integer  iat,iap,izt,izp, a, z,n, target jieuts,projjieuts,ie,ne,i 

PI-4.*DATAN(1.D0) 
target_neuts=iat-izt 
proj_neuts=iap-izp 
nucleon_mass=.93  86 
r0=1.26 

c  proton_factor  will  be  used  to  convert  neutron  cross  sections  to  protons 
proton_factor=(real(izt+izp))/(real(target_neuts+proj_neuts)) 

c  sigmaO  is  a  parameterized  absorption  cross  section  in  two  forms: 
c  the  nucleon-nucleus  form  is  due  to  Wilson  et.  al.  1988 
c  the  nucleus-nucleus  form  is  due  to  Townsend  et.  al.  1986 
if((iat.eq.  1  ).or.(iap.eq.  1))  then 
f=l-0.62*Exp(-tnnlab/200)*Sin(10.9*tnnlab**(-.28)) 
if(iat.eq.l)  then 

sigma0=f*45.*real(iap)**.7*(l. +0.01 6*Sin(5. 3-2. 63*Log 
&(real(iap)))) 
else 

sigma0=fl'45.*real(iat)**.7*(l.+0.016*Sin(5.3-2.63*Log 

&(real(iat)))) 


end  if 
else 

beta=l.+5./tnnlab 

delta=0.2+l./real(iap)+l/real(iat)-0.292*Exp(-tnnlab/ 

&792.)*Cos(0.229*tnnlab**0.453) 

inner=(real(iap)*  *  ( 1  ,/3  .)+real(iat)*  *(  1  ,/3  .)-delta) 
sigma0=10.*pi*(r0**2)*beta*(inner**2) 
end  if 

c  now  start  looping  through  energies  and  values 
do  10  ie=l,ne 
elab=.001*eel(ie) 
p_total=linvt(ie) 
print  ' 

print  *,'ie',ie,'p_total',p_total 

c  now  set  up  cases  for  each  fragment  of  interest 
c  l=deuteron,  2=triton,  3=helion,  4=alpha 
c  pO  is  the  coalescence  radius  for  the  resultant  nuclei  (MeV/c) 
c  these  are  approximations  of  those  used  in  Gutbrod  et.  al 
c  deut_p0=125 
c  tritjp0=129 

c  helion_p0=129  (observed  production  of  helions  is  It  tritons) 
c  alpha_p0=140 

c  these  values  are  from  Awes  et.  al.  (1981);  helion  interpolated 
c  deutjp0=170 
c  trit_p0=2 1 5 
c  helionjp0=240  (approx) 
c  alpha_p0=270 

c  Nagamiya  et.  al.  (1981)  suggested  pO  of  ~90 

do  20  i=l,4 
select  case(i) 
case(l) 
z-1 
n=l 
a=2 
p0=.09 
case(2) 
z=l 
n=2 
a=3 
p0=.09 
case(3) 
z=2 
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11=1 
a=3 
p0=.09 
case(4) 
z=2 
n=2 
a=4 
p0=.09 
end  select 

gamma=l  +elab/nucleon_mass 

k(i)=(a/(z*n))*proton_factor**z*(4.*pi*p0**3./(3.*sigma0* 
&nucleon_mass))  *  *  (a- 1 ) 

if(ie.eq.l)print  *,'i',i,'k(i)',k(i) 

20  continue 

c  now  set  up  cases  for  each  fragment  of  interest 
c  r=deuteron,  s=triton,  t=helion,  u=alpha 

r=k(l) 

s=k(2) 

t=k(3) 

u=k(4) 

q=s+t 


c  this  section  solves  the  following  cubic  equation: 
c  p_total=p  obs  1  +r*p_obs*  *2+q*p_obs*  *3  for  p_obs  given  p_total 
c  neglecting  the  (very  small)  alpha  term.  Solving  the  fourth  order 
c  polynomial  is  not  unreasonably  complex  in  itself,  but  fortran  doesn't 
c  like  to  work  with  the  cube  root  of  a  negative  number  -  wants  to  make  a 
c  complex  number  out  of  it.  The  error  is  very,  very  small  until  you  get 
c  to  extremely  large  (-10**7  mb)  cross  section,  then  grows  slowly 
parta=qsqrt(27.*q**2*p_total**2-r**2*(l.+4.*r*p_total)+q* 
&(4.+18.*r*p_total)) 

partb=(2.*r**3-27.*q**2*p_total+3.*q*(-3.*r+sqrt(3.)* 

&parta))**(l./3.) 

p_obsl=-(2.*r-(2.*2.**(l./3.)*(3.*q-r**2))/partb+2.**(2./3. 

&)*partb)/(6.*q) 

c  these  are  the  two  reduction  factors  -  rfl  is  based  on  the  calculated 
c  coefficients  using  the  total  production  cross  sections,  rf2  is  based 
c  on  the  parameterized  coefficients  used  to  calculate  the  observed  cross 
c  sections.  At  present,  rf2  is  used.  It  is  possible  to  use  the 
c  calculated  coefficients  in  the  solution  for  the  observed  cross  sections, 
c  However,  the  results  should  differ  only  slightly, 
rfl  =(  1  ,-r*p_total-(s+t)*p_total*  *2-u*p_total*  *  3) 
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rf2=p_obs  1  /p_total 
print  *,'rf2',rf2 
ds_total(ie)=rf2  *  ds_total(ie) 
print  *,'p_obsl',p_obsl,'d  after  (rf2)',ds_total(ie) 
10  continue 


return 

end 


Appendix  C 


Additional  Plots 


75  MeV  per  nucleon  C+Al 


Table  C.l  Chi-squared  goodness  of  fit  check  for  previous  and  current  models.  Bold  type 
indicates  best  fit. 
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Figure  C.l  Comparison  of  results  for  75  MeV  per  nucleon  C+Al  at  30°  lab  angle 
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Figure  C.2  Comparison  of  results  for  75  MeV  per  nucleon  C+Al  at  45°  lab  angle 
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Figure  C.3  Comparison  of  results  for  75  MeV  per  nucleon  C+Al  at  90°  lab  angle 
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